Single-particle spectral function of the Holstein-Hubbard bipolaron 
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The one-electron spectral function of the Holstein-Hubbard bipolaron in one dimension is stud- 
ied using cluster perturbation theory together with the Lanczos method. In contrast to other 
approaches, this allows one to calculate the spectrum at continuous wavevectors and thereby to 
investigate, for the first time, the dispersion and the spectral weight of quasiparticle features. The 
formation of polarons and bipolarons, and their manifestation in the spectral properties of the 
system, is studied for the cases of intermediate and large phonon frequencies, with and without 
Coulomb repulsion. A good agreement is found with the most accurate calculations of the bipo- 
laron band dispersion available. Pronounced deviations of the bipolaron band structure from a 
simple tight-binding band are observed, which can be attributed to next-nearest-neighbor hopping 
processes. 

PACS numbers: 63.20.Kr, 71.27.+a, 71.38.-k, 71.38.Mx 
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I. INTRODUCTION 

In recent years, angle-resolved photoemission spec- 
troscopy (ARPES) has proved to be very helpful in 
obtaining information about the electronic states of 
strongly correlated systems. While a lot of data is avail- 
able from experiments, reliable theoretical calculations of 
the one-particle spectral function — which can often be re- 
garded as being proportional to the ARPES spectrum — 
for popular models of, e.g., the Hubbard or t-J type, 
are usually very demanding. As a consequence, many 
interesting problems of condensed matter physics have 
not been investigated systematically regarding their spec- 
tral properties in a satisfactory way. Among them is the 
"bipolaron problem" of two electrons, which can form a 
bound state even in the presence of strong Coulomb re- 
pulsion if they are coupled to phonons. Despite the long 
history of this problem, bipolaron formation is still sub- 
ject of ongoing discussion due to its potential role, e.g., 
in high-temperature superconductors^ and manganites, 2 
two classes of materials studied extensively over the past 
decade. 

Most existing results for the Holstein-Hubbard (HH) 
model considered here have been obtained using exact 
diagonalization (ED). Apart from a systematic error due 
to the necessary truncation of the Hilbert space, this 
method gives exact results for the one-electron spectrum, 
but is restricted to rather small systems, especially for 
small phonon frequencies and/or strong electron- phonon 
coupling. Consequently, it is difficult to study the dis- 
persion of the spectral peaks throughout the Brillouin 
zone. To overcome this limitation, we employ here clus- 
ter perturbation theory (CPT), extending the recent ap- 
plication to the Holstein model with one electron^ In 
contrast to ED, CPT permits to calculate the spectral 
function for continuous wavevectors. Moreover, finite- 
size effects are strongly reduced compared to direct diag- 
onalization of small clusters, and results are much more 
realistic than previous work based on, e.g., a two-site 
systemA^ CPT becomes exact in the weak- and strong- 



coupling regime, and has been successfully applied also to 
other problems iSi2iiSiiiii2iiii4 A review of cluster meth- 
ods for strongly correlated systems is by Maier et al. 15 
Here we would merely like to point out the recent ap- 
plication of the dynamical cluster approximation to the 
half-filled Holstein model by Hague^ 

In this work, we study in detail the formation of po- 
larons and bipolarons, and its dependence on phonon 
frequency, electron-phonon and electron-electron inter- 
action. To test the reliability of our results, and to in- 
terpret the quasiparticle (QP) features in the spectra, 
comparison is made to ED, as well as to the most ac- 
curate approach to the one-dimensional HH bipolaron 
currently available, namely the variational diagonaliza- 
tion method ml Moreover, we also investigate the form of 
the bipolaron band dispersion and compare it to that of 
models with nearest-neighbor- and next-nearest-neighbor 
hopping. 

This paper is organized as follows. We begin with a 
review of the HH bipolaron in Sec. [n] I n Sec. II I II we 
discuss some details of the application of CPT. Results 
are presented in Sec. II VI and, finally, Sec. IVl contains our 
conclusions. 



II. THE HOLSTEIN-HUBBARD MODEL 



The HH model is defined by the Hamiltonian 



H 



E 



gYMbt+b^ + U 



(1) 



Here c\ a (c icr ) and b\ (b t ) create (annihilate) an electron 
of spin cr, and a phonon of energy u> (h = 1) at site i, re- 
spectively, and Hi = J2 a n ia with n ia = Ci a c ia- The first 
two terms correspond to the kinetic energy of the elec- 
trons and the kinetic and elastic energy of the phonons, 



2 



respectively. The electron-phonon (el-ph) and electron- 
electron (el-el) interaction is described by the third and 
forth term. We have three model parameters, namely 
the amplitude for nearest-neighbor hopping t, the phonon 
frequency u>, the el-ph coupling constant g, and the el-el 
interaction strength U > 0. For {7 = 0, Eq. JJJ is identi- 
cal to the pure Holstein model, while for g — we recover 
the Hubbard model. We introduce the commonly used 
dimensionless coupling constant A = 2g 2 /(ujW), where 
W = 4tD is the bare bandwidth in D dimensions. We 
further define the dimensionless parameters U = u>/t and 
U = U/t, and express all energies in units of t. Conse- 
quently, the independent parameters of the model are uj, 
U, and A. In the sequel, we shall also use the polaron 
binding energy Ep — XW/2, which emerges as a nat- 
ural parameter from the Lang-Firsov transformation. 18 
Finally, the lattice constant is taken to be unity. 

Owing to the complexity of even the two-electron prob- 
lem, we will not discuss effects of bipolaron-bipolaron 
interaction here. An investigation of the latter, which 
will definitely play an important role in real materials, 
requires a study of the HH model with many electrons 
(see, e.g., Ref.ll3and references therein), which is beyond 
the scope of our method in its present form. 

There exists a considerable amount of work on the HH 
bipolaron, although it is by far not as well understood 
as the simpler one-electron case. In the sequel, we re- 
strict our discussion to new developments in the field. A 
very complete review of earlier work has been given by 
Alexandrov and MottiSS 

While the pairing of electrons in momentum space can 
be accurately described by Migdal-Eliashberg theory 21 
for weak-enough coupling, no reliable theory is avail- 
able for the formation of bipolarons — corresponding to 
pairing of electrons in real space — at intermediate to 
strong el-ph interaction. In recent years, progress was 
made using either variational approache o 22 ' 23 i 24 ' 25 i 26 ' 27 i 28 
or, more importantly, unbiased numerical studies based 
on ED^ 4 ! 5 ! 6 ! 29 ! 30 variational diagonalizationt 17 ! 31 ! 32 ! 33 
DMRG, 34 and QMC, 35 i 36 i 37 The ED and DMRG calcu- 
lations were restricted to rather small systems consist- 
ing of two; 4 i 5 i 6 four^ six ml eight 30 ! 31 or twelve sites ^ 
while the methods of Refs . UtI l3a andl36l are almost free of 
finite-size effects. The larger number of phonon states re- 
quired to obtain converged results makes numerical stud- 
ies with ED methods even more challenging than for a 
single electron, especially for small phonon frequencies. 

Since the HH model represents a simplified description 
of the situation in real materials, it is highly desirable 
to study more complex models. To this end, it is in- 
teresting to note that the QMC methods of de Raedt 
and Lagendijk 3 ^ and Macridin et alm& may be general- 
ized to include dispersive phonons. Furthermore, both 
approaches can be applied to models with long-range 
Coulomb interaction; 3 ^* 3 ^ similar to the work of Bonca 
and Trugman^ Finally, bipolaron formation in a model 
with Jahn- Teller modes — as present, e.g., in perovskite 
manganites — has recently been investigated by Shawish 



et alM 

To discuss bipolaron formation in the HH model, we 
have to distinguish between two cases. The two electrons 
can either have the same or opposite spin, which leads to 
a singlet or triplet state, respectively. We consider these 
possibilities separately. 



1. Singlet state 

For two electrons in a singlet state, the formation of a 
bound bipolaron state in the absence of Coulomb inter- 
action originates from the fact that the potential well — 
arising from a displacement of the oscillators — around an 
occupied lattice site deepens in the presence of a second 
electron. This may easily be seen in the atomic limit 
t = 0, using the Lang-Firsov transformation^ On differ- 
ent lattice sites, each electron gains an energy — Ep by 
distorting the lattice, whereas the energy shift becomes 
—4Ep if both particles occupy the same site (small or on- 
site bipolaron). For t ^ 0, the competition between the 
kinetic energy of the electrons on the one hand and the 
displacement or lattice energy on the other hand deter- 
mines the cross over from a state with two weakly bound 
polarons, sometimes also referred to as a large bipolaron, 
for A < A c to a small bipolaron for A > A c , where A c de- 
notes the critical value of the el-ph coupling. In Sec. ITT1 
A has been defined as A = 2Ep/W, i.e., as the ratio of 
the energy gain due to polaron formation to the kinetic 
energy of a free electron. While A c = 1 in the adiabatic 
regime for the small-polaron cross over in the model with 
one electron (see, e.g., discussion in Ref. l39j) . here we ex- 
pect A c = 0.5 (for uj <c 1) due to the energy gain of —2Ep 
per electron compared to — Ep in the single polaron prob- 
lem. This is well confirmed by the calculations of Wellein 
et al.^SL who find a strong decrease of the kinetic energy 
near A = 0.5 for uj = 0.4. 

For uj 3> 1, the lattice energy becomes important 
since the trapping of the carriers requires a sizable lat- 
tice distortion. This gives rise to the additional condition 
2 y/ Ep I uj > 1 for a small bipolaron 4°« Similar to the one- 
electron problem, the cross over is very gradual in the 
nonadiabatic regime«2°» The correlation or binding of the 
two electrons depends crucially on the phonon frequency, 
since the latter determines the maximum distance across 
which the two particles feel an attractive interaction due 
to the phonons. Up to second order in g, this coupling is 
given by 

U cS (e) = g 2 V ph (q,e) = -^^ 1 , (2) 

where Z? p h(q,e) denotes the phonon propagator. Equa- 
tion J2Jl reveals that the energy-dependent interaction is 
attractive for e < uj, and becomes instantaneous in the 
antiadiabatic limit ui — > oo where U c s = —2g 2 /uj. Hence, 
the binding always decreases with increasing phonon 
frequency^ 
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For U > 0, there is a competition between the re- 
tarded, attractive interaction mediated by the phonons 
and the instantaneous, repulsive Hubbard interaction. 
Consequently, a state with two unbound polarons — 
stabilized by the onsite repulsion — can exist for suffi- 
ciently weak el-ph coupling. 30 This is in contrast to the 
extended HH model with long-range interaction, in which 
a bipolaron state is formed irrespective of the value of 
t/4i The effective el-el interaction in the HH model de- 
termining the nature of the bipolaron state is 



U eS = U- 2E P . (3) 



From this result, which can be obtained either from the 
generalization of Eq. J2Jl to U ^ in the limit u> — > oo, 
or in the antiadiabatic strong-coupling limit ^ one may 
be tempted to expect a bipolaron state to exist only 
for U c {f < 0, i.e., if there is a net attractive interaction 
between the particles. While this is true for the effec- 
tive Hubbard model onto which the HH model maps in 
the antiadiabatic strong-coupling limit, a consideration 
of virtual hopping processes leads to the less stringent 
condition U < 4Ep H The energy gain due to virtual ex- 
change processes of two electrons on neighboring lattice 
sites — not suppressed by strong el-ph coupling — permits 
the formation of a weakly bound intersite bipolaron with 
the two electrons most likely to reside on neighboring 
lattice sitesiii2& A phase diagram for bipolaron forma- 
tion as a function of A and u; in one dimension has been 
presented by WeiBe et al£& Eventually, for sufficiently 
strong el-ph coupling 2Ep > U, the effective onsite po- 
tential U e f{ becomes attractive, and a small bipolaron is 
formed. 

Starting from a small bipolaron, a cross over to an in- 
tersite bipolaron takes place when the Coulomb interac- 
tion becomes large enoughii 7 * 3 ^^^ The intersite bipolaron 
has a much smaller effective mass than the small bipo- 
laron and may therefore also exist as a mobile carrier in 
real materials^ In the adiabatic limit to = 0, the onsite- 
intersite bipolaron transition has been shown to be of 
first order but for finite phonon frequencies it is ex- 
pected to happen in a more gradual way because of re- 
tardation effects, in agreement with recent calculations. 17 
Estimates for the region of existence of the intersite bipo- 
laron state for Zo = 1 are U < 2Ep for weak coupling, and 
U < 4Ep for strong el-ph coupling, 1 and phase diagrams 
in the (U, A)-plane have been reported in onei 7 - and two 
dimensions. 36 While the above conditions are quite accu- 
rate in the nonadiabatic regime W > 1, the case 57 <C 1 
remains an open problem. 

Finally, the physically most interesting regime, which 
is unfortunately also the most difficult case to treat the- 
oretically, is defined by To -C 1 , and a Coulomb repulsion 
at least as large as the attractive interaction due to the 
el-ph coupling. 



2. Triplet state 

For two electrons of the same spin, the Pauli principle 
forbids double occupation of a site. In principle, a bound 
state may be formed with the two particles being located 
on different lattice sites. While two electrons can lower 
their energy by sharing a lattice distortion in the large 
bipolaron regime, especially for small phonon frequen- 
cies, the exchange process stabilizing the singlet intersite 
bipolaron state at intermediate-to-strong el-ph coupling 
and U > is not strong enough to bind two polarons in 
a triplet intersite stated Furthermore, for U < oo, the 
ground-state energy of the triplet state is always larger 
than for a singlet state because two particles with par- 
allel spin cannot occupy the same k = energy level. 
Finally, the singlet and triplet states become degenerate 
in the limit U — > oo. 



III. METHOD 

As mentioned above, here we use CPT in combination 
with the Lanczos recursion method^ Details about the 
application to el-ph problems have been given in Ref.H 
henceforth also referred to as I. The major difficulty we 
are facing in the present case is the larger number of 
phonon states needed to obtain converged results. From 
a physical point of view, this is not surprising since each 
of the two electrons will create a lattice distortion or 
phonon cloud, whereas there is only one dressed parti- 
cle (polaron) in the case of the Holstcin polaron consid- 
ered in I. However, in addition to the simple doubling of 
the number of particles, it has been shown by previous 
authors22i224^4 that multiphonon states play a more im- 
portant role for the bipolaron as a result of the phonon- 
mediated binding. 

ED (and also CPT) for el-ph systems is affected both 
by finite-size effects and the truncation error due to the 
restricted number of phonon states kept in calculations. 
Obviously, if one used very small clusters, good conver- 
gence with respect to the phonons could be achieved even 
for strong el-ph coupling. On the other hand, for small 
numbers of phonon states, rather large clusters can be 
studied. The approach which has been widely used in 
the past is to require the truncation error, e.g., of the 
ground-state energy, to be smaller than a certain limit, 
and to use the maximal cluster size which can be handled 
for this number of phonons. Here, an additional challenge 
arises form the fact that the diagonalization of the cluster 
has to be performed for open boundary conditions. 3 Con- 
sequently, one cannot exploit translational symmetry to 
reduce the dimensionality of the Hilbert space. However, 
this drawback is clearly outweighed by the advantages of 
CPT outlined below. 

We would like to briefly discuss some interesting fea- 
tures of CPT. The method is based on a breakup of the 
infinite lattice into clusters of N sites, say£ The one- 
electron cluster Green function, denoted here as G a b,a 
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[see Eq. (4) in I] , of the model under consideration is cal- 
culated for one of these (identical) clusters using open 
boundary conditions. This can be done, e.g., using ED 
or analytical approaches £ The hopping between adjacent 
clusters is then treated as a perturbation to obtain the 
Green function of the original system, Go-(k, e)£ A basic 
limitation of the theory is that the Hamiltonian must not 
contain any nonlocal interactions, except for one-electron 
terms. Additionally, CPT in its present form can only be 
used to calculate one-particle Green functions^ There- 
fore, interesting observables such as, e.g., el-el correlation 
functions or transport properties are not yet available. 

From the nature of the approximation made, it is clear 
that CPT will work particularly well if the local inter- 
actions dominate the physics of the system, i.e., for the 
case of the HH model (Q, if g, U ^> t. This point will 
be illustrated in Sec. IIVI The method becomes exact 
in the atomic limit t = 0, for nonintcracting electrons 
(g,U = 0), as well as for N = ooA The quality of the 
results obtained with CPT has been tested for several 
modelsj£*2i2ii2iiiii2ii&ii and a very good agreement with 
other work has been found. In I, we pointed out the oc- 
currence of finite-size effects which show up as additional 
peaks in the corresponding one-electron spectral func- 
tion. The weight of the latter reduces quickly as N —> oo, 
so that the spectrum is not affected significantly. 

One of the most important advantages of CPT is that it 
gives, in principle, results for an infinite system, although 
the approximate treatment of the intercluster hopping 
introduces some finite-size effects, which can be system- 
atically reduced by increasing the cluster size. As a con- 
sequence, the one-electron Green function can be calcu- 
lated for any wavevector in the Brillouin zone, even for 
N = 1. This allows one to study the dispersion of the 
QP peaks, in strong contrast to standard ED methods on 
finite clusters, for which only a few points in momentum 
space are accessible, owing to the rather small values of 
N w 2-20 usually used. 

Since we only consider the one-dimensional HH model 
in the sequel, we shall adopt the notation accordingly. 
We are interested in the one-electron Green function 

where ||) denotes the ground state with one electron 
of spin down, and a =t,|- Equation Q only contains 
the inverse photoemission part of the total one-electron 
Green function. In the case of G-\, the second part — 
corresponding to photoemission (PE) — vanishes, since 
there is no f electron in the ground state. Moreover, 
for Gi, PE is identical to the spectrum of a single po- 
laron, which has been studied in detail in I. The situa- 
tion would be different if we started with a two-electron 
(singlet or triplet) ground state. Then, the PE part of 
the one-electron spectral function also contains valuable 
information. However, due to limited computer memory, 
such computations involving three-electron states are not 
possible with the code used here. 



In Eq. Q, we have omitted the energy Eq of the 
ground state J.}, which usually enters in the form H—Eq, 
to permit direct comparison with the singlet bipolaron 
band dispersion E^(k) in Sec. IIVI The one-electron 
spectral function is related to the Green function 10} via 

A a (k, e) = -7T -1 lim lmG a (k, e + irj) . (5) 

To calculate the cluster Green function by ED, a trun- 
cation of the phonon Hilbert space is necessary, and we 
use the same truncation scheme as in I. The number of 
phonon states iV p h will be chosen so as to push the trun- 
cation error A = |£ n (AT ph + 1) - E^ 1 (N ph )\/\El l (N ph )\ 
of the energy Eq of the two-electron ground state If J.) 
below 10~ 4 . The use of to monitor convergence with 
respect to 7V p h comes from the observation that — for the 
same number of phonons — the truncation error of the 
latter is always smaller than for the triplet state \[[). 
This may be ascribed to the fact that for two electrons 
of the same spin, no bound onsite or intersite state ex- 
ists. In particular, there will be no large local lattice 
distortion surrounding an onsite bipolaron, the descrip- 
tion of which requires a significant number of phonons. 
In previous work on the HH bipolaron, using ED with pe- 
riodic boundary conditions^*Si22i2£ the truncation error 
was usually smaller than 10 -6 . However, these meth- 
ods were restricted to only a few k vectors. Furthermore, 
our calculations show that even a relative error A = I CP 4 
ensures satisfactory convergence of the one-electron spec- 
trum. The smaller number of phonon states enables us 
to use larger clusters and thereby significantly diminish 
finite-size effects since, within the CPT, even an increase 
N — > N + 1 noticeably improves the results. Once the 
cluster size has been fixed, we use the maximal possi- 
ble number of phonons. The accuracy A varies for the 
different calculations and will be reported in each figure. 

In its present form, our method is restricted to the 
nonadiabatic regime ZU > 1, except for weak el-ph cou- 
pling. To study smaller phonon frequencies — relevant to, 
e.g., to transition metal oxides — a combination with vari- 
ational diagonalization techniques, or the use of shared- 
memory systems would be necessary. As in I, we restrict 
our calculations to the spectral function, which is the 
most fundamental quantity that can be obtained from 
CPT£ 



IV. RESULTS 

The one-electron spectral function of the problem con- 
sidered here has been calculated before using ED 4 *^ and 
DMRG, 34 both in one dimension. However, results were 
only given for k = 0, and for very small systems with 
N = 2 and N — 6, respectively. With the above meth- 
ods, and for periodic boundary conditions, the spectral 
function (JjJJ can be evaluated for N different wavevec- 
tors, of which only N/2 + 1 are physically nonequiva- 
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lent. This makes it difficult or even impossible to study 
the dispersion of QP features. Recently, a parallelized 
DMRG code has been developed,*^ which allows studies 
of one-dimensional Holstcin models on very large clusters 
even at half filling^ However, the calculation of spectral 
functions within DMRG is very time-consuming, since it 
has to be done separately for each point on the energy 
axis. Several authors have also calculated dressed spec- 
tral functionsj£*§i2i with the fermion operators in Eq. (0J 
replaced by their Lang-Firsov transformed (i.e., dressed) 
counterparts, as well as pair spectral functions^ The 
corresponding spectra show a simplified structure in cer- 
tain regimes, indicating that polarons and bipolarons are 
"good" QP's for these parameters. 

De Mello and Ranninger— have pointed out that to 
study the cross over between polarons and bipolarons 
it is, in general, necessary to investigate both, photoe- 
mission and inverse photoemission. This can easily be 
understood by considering electron emission from the 
two-electron singlet ground state, i.e., the Green function 
(til c I|( e — H)~ lc ki ITD- Depending on the parameters, 
|t|) may either consist of two weakly bound polarons or 
a bipolaron. Consequently, photoemission spectra will 
only show a single QP band. In contrast, the Green func- 
tion with a =| corresponds to adding an \ electron 
to the one-electron ground state ||). For example, the 
additional particle can either go into the ground state to 
form a bipolaron, or into an excited polaron state. In 
general, we therefore expect two QP bands in the spec- 
tral function, whose weights, positions and widths vary 
with ZD, U and A. 

As we will compare our findings with the variational 
diagonalization method (VDM) of Bonca et nl.m^M we 
would like to comment on the accuracy of the latter. The 
problem is defined on an infinite system, so that the ap- 
proach is free of boundary finite-size effects. However, the 
method involves a variationally determined Hilbert space 
with two variational parameters, namely the maximal al- 
lowed distance between electrons and phonons, and be- 
tween the two electrons, respectively. For the bipolaron 
problem under consideration, the limiting parameter in 
the regime ZD > 1 is the maximum distance iVh between 
the two electrons. The results presented here have been 
obtained using < 18. While the method gives very ac- 
curate results — with errors smaller than the linewidth in 
the figures — for the case of a small bipolaron (U <C 2i?p), 
it is less reliable (relative errors < 1%) for strong on- 
site repulsion U 3> 2Ep favoring two weakly bound po- 
larons, similar to ED and CPT. Due to additional towers 
of phonon excitations that are located in the neighbor- 
hood of the electron sites, the method achieves good con- 
vergence in the small bipolaron regime even for strong 
coupling. Nevertheless, the adiabatic regime ZU <C 1 rep- 
resents a difficult problem, as is the case for other ap- 
proaches. Finally, as in CPT, results may be obtained at 
any wavevector. 

We shall see below that there is a close correspon- 
dence of the QP bands in the spectra to the polaron and 



bipolaron dispersion relations denoted here as E^{k) and 
E^(k), respectively. The notation (fc) is convenient, 
but there is no spin-dependence in the case of a single 
electron, i.e., £ T (fc) = E L (k). Results for E n (k) have 
been reported by Wellein et alm& and Weifie et al^J^ 
However, in contrast to A a (k,e), E^(k) and E^(k) do 
not reveal the spectral weight of the corresponding QP's. 
Nevertheless, the comparison with the spectra will yield 
valuable insight and serve as a test of the CPT results. 
Moreover, a direct calculation of energy bands does not 
suffer from the restricted energy resolution of CPT due 
to the use of a smearing parameter [Eq. ||SJ|]. 

Owing to the limitations regarding the number of 
phonon states, we shall only show results for ZD > 1. To 
be more specific, we consider two values of the adiabatic 
ratio, namely ZD = 4 and ZD = 1. For ZD = 4, the spectra 
will turn out to be relatively simple, and we are able to 
study even strong el-ph coupling. Consequently, we start 
with a discussion of the antiadiabatic regime, and then 
move on to the more difficult case ZU = 1. 



A. Antiadiabatic regime 

In this section, we restrict the discussion to U = 0, 
while the influence of Coulomb repulsion will be studied 
below. Figure ^ shows the evolution of the one-particle 
spectrum with increasing el-ph coupling. Here and in 
subsequent figures, solid lines represent results for 
and dashed lines correspond to Aj_. 

For [7 = 0, two electrons of opposite spin always form a 
bipolaron state for any A > 0. At weak coupling A = 0.5 
[Fig.^a)], j4| exhibits two well visible bands, as well as 
an incoherent part centered at e ~ 0. To understand the 
nature of the coherent excitations, we have also included 
in Fig.^the bipolaron band dispersion E'^(k) (solid ver- 
tical line), calculated by the VDM^i The latter fits well 
the low-energy band, with the minor deviations at inter- 
mediate k — i.e., the splitting of the low-energy peak into 
several small satellites — being finite-size effects, as has 
been verified by calculations on smaller and larger clus- 
ters for a smaller number of phonon states (not shown). 
A more detailed discussion of finite-size effects will be 
given below for ZD = 1 . 

Even at weak coupling A = 0.5, the bipolaron band 
already has a relatively small width of W /W « 0.37 
compared to the free-electron value W . Moreover, the 
spectral weight of the lowest-energy peak, obtained by in- 
tegration over the CPT spectrum, decreases significantly 
from about 0.68 at k = to about 0.08 at k — it. At the 
same time, the weight contained in the incoherent part of 
the spectrum increases with increasing k. This behavior 
is very similar to the single-electron case£ 

We now turn our attention to the second, higher-lying 
band appearing in Fig. ^a). From the general discus- 
sion in Sec. [HI we expect that it corresponds to an ex- 
cited state with two polarons. We therefore compare it 
to the energy of two independent polarons in an infinite 
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FIG. 1: (color online) Spectral functions Af(k,e) (solid lines) and Aj(fc, e) (dashed lines) calculated with CPT for different 
values of the el-ph coupling A, using r\ = 0.05i [see Eq. 0]. All other parameters as indicated in the figures. The truncation 
errors are A < (a) 5.3 x 1CT 6 , (b) 1.0 x 1CT 5 , (c) 2.4 x 1CT 7 , (d) 6.7 x 10~ 6 (see text). The vertical lines correspond to VDM 
results for the polaron and bipolaron band dispersions -B T (0) + E'(k) (dashed) and E^ l {k) (solid), respectively^ 



system. Since Aj describes the process of adding an elec- 
tron with momentum k to the one-polaron ground state 
with energy E^(0) [cf. Eq. Q], we show in Fig. QJa) 
the band dispersion E^(0) + E^ (k) (dashed vertical line). 
The comparison with the spectral function yields a very 
good agreement at intermediate and large k, while there 
are some discrepancies at small momenta. A density plot 
of Af (Fig. EJ reveals more clearly that the two coher- 
ent bands hybridize and repel each other near the point 
where they would be degenerate, giving rise to an upper 
band with inverscd dispersion at small k. The situa- 
tion is similar to the hybridization of the coherent and 
incoherent parts in the one-electron case occurring for 
\E^(k) — E^(0)\ ~ ZJ (see I). Of course, such effects are 
absent in the band dispersion of a system with two inde- 
pendent polarons. Since the residual interaction between 
two polarons vanishes in the limit N — > oo, the hybridiza- 
tion visible in the CPT spectrum may be attributed to 
finite-size effects. The latter originate from the fact that 
within CPT, translational symmetry is broken by treat- 
ing inter- and intracluster hopping differently, and only 
approximately restored afterward. 

The spectral function A±, also shown in Fig.^a), con- 



tains a coherent band at low energies, and an incoherent 
part which is very similar to that of Af . Well away from 
k = 0, the coherent peaks in A± follow closely the polaron 
band in A-p . Thus the excited two-polaron state of the 
system with two electrons of opposite spin is very similar 
to the ground state of the system with two electrons of 
the same spin. Near k = 0, the spectral weight of the 
low-energy peak in A| is small (~ 0.08) compared to the 
polaron peak in Aj (« 0.2). This is a result of the fact 
that two polarons with the same spin cannot occupy the 
same k = state. The picture changes at larger mo- 
menta, where both bands have similar weight, although 
the sharp peaks in A^ are higher than the broadened 
features in A^. 

With increasing el-ph coupling, the low-energy bipo- 
laron band becomes even narrower until it is virtually flat 
at A = 1.5 [Fig.njc)]. Here, the two conditions for a small 
bipolaron (Sec.|nJ are identical to A > 0.5. Consequently, 
finite-size effects are very small in Figs.^b) - (d), as con- 
firmed by the excellent agreement of the CPT data with 
the results for E^(k). The reduction in bandwidth with 
increasing A is accompanied by a loss of spectral weight. 
For k = 0, the latter decreases from the value 0.68 at 
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FIG. 2: (color online) Density plot of the spectral func- 
tion Af(k,e) for ZJ = 4.0, U = 0, and A = 0.5, as shown 
in Fig. Ua). The symbols correspond to VDM results for 
E^ (0) + E^ (k) (squares) and E^(k) (crosses), respectively, 41 * 



A = 0.5 given above to about 0.10 at A = 2.0. Both the 
narrowing and the loss of weight indicate a significant 
increase of the effective bipolaron mass. 

While the polaron band lies relatively close to the bipo- 
laron band at A = 0.5 [Fig. da)], the increase of the 
coupling leads to a clear separation, and to a downward 
shift of both bands proportional to the polaron binding 
energy Ep. In the antiadiabatic strong-coupling regime 
of Fig.[TJd), the energy gap between the two bands is well 
described by the atomic- limit value 2Ep = 8t. Similar to 
A = 0.5, the two-polaron band dispersion E^(0) + E^{k) 
agrees well with the polaron band in the spectra, with 
some differences being visible near k = 0. Interestingly, 
in Fig.^c), there is a mixing of the bipolaron state with 
one phonon excited, which lies an energy to = At above 
the lowest band, and the two-polaron excitation. 

The polaron band also narrows with increasing el-ph 
coupling (see also I). However, the effect is much smaller 
than for the bipolaron band. Additionally, the spectral 
weight of the k — polaron peak in increases from 
about 0.20 at A = 0.5 to about 0.32 at A = 2.0. This 
may be explained by the fact that for weak coupling 
[Fig. ma)], some of the weight of the polaron state is 
contained in the large low-energy feature. Calculations 
for a single electron and the same parameters show that 
the spectral weight of the polaron at k = decreases 
from about 0.86 (A = 0.5) to about 0.52 (A = 2.0). Since 
the spectral weight is, to a very good approximation, 
equal to the inverse of the effective mass of the Holstein 
polaron^ these results indicate that the polaron mass 
does not increase at the same rate as the bipolaron mass 
with increasing coupling, as reflected by the correspond- 



ing changes in bandwidth in Fig. m Finally, we also find 
a comparable reduction of spectral weight for the two- 
polaron band in A± from about 0.08 (A = 0.5) to about 
0.04 (A = 2.0) at k = 0. 

To conclude the discussion of the case ZJ = 4, we would 
like to underline the enormous advantage of CPT in the 
strong-coupling regime. It permits us to perform calcu- 
lations on a very small cluster (N = 4) — sufficient to ob- 
tain well-converged results — but still yields the spectral 
function at any desired k. 



B. Intermediate phonon frequency 

In the preceding section, we have investigated in detail 
the signatures of polaron and bipolaron states in the one- 
particle spectrum for ZJ = 4. Owing to the large energy 
of phonon excitations, most of the spectral weight resides 
in the corresponding bands, allowing a fairly easy iden- 
tification. We now consider the case ZJ = 1 , which turns 
out to be more difficult to study numerically and to in- 
terpret. Nevertheless, work in the regime ZJ < 1 is highly 
desirable to understand many interesting strongly corre- 
lated systems such as, e.g., the manganites. Although 
the latter are usually characterized by u> <C t, quantum 
effects are already visible for uj = t. As a consequence, 
previous author ^J^^^i^ have often focused on this 
case, which is numerically much easier to tackle than the 
region ZJ <C 1. The case ZJ <C 1 has been considered, e.g., 
by Wellein et aL 30 and Weifie et nlMJL While the dis- 
cussion for uj = 4 was restricted to U = 0, here we shall 
also take into account a finite Coulomb repulsion. 



1. (7 = 

Since converged results for ZJ = 1 require more phonon 
states than for ZJ = 4, we have slightly reduced the cluster 
sizes in our calculations. Consequently, finite-size effects 
are larger, as discussed below. Moreover, we are not able 
to reach the strong-coupling regime but instead restrict 
the range of A to 0.5-1.25. 

Figure 13 contains the one-particle spectra for U = 0. 
In principle, for A = 0.5, the results look quite simi- 
lar to Fig. 01a). However, the spectral weight of the 
two coherent bands is much smaller, as a consequence 
of the increased importance of incoherent excitations for 
ZJ = 1. In particular, the weight of the latter is strongly 
enhanced at large k, so that the bands are no longer easy 
to identify. Therefore, and because of the strong mixing 
of the bands with coherent and incoherent excitations, 
it becomes difficult to accurately determine the spectral 
weight by integration over the CPT spectra. 

We see from Fig. [3] that the bipolaron bandwidth is 
much smaller for ZJ = 1 (W /W « 0.1) than for ZJ = 4 
[Fig. ma)], despite the fact that the value of A is the 
same in both cases. Hence, the effect of el-ph interaction 
on the bipolaron mass is much more pronounced in or 
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FIG. 3: (color online) Spectral functions Af(k,e) (solid lines) and Ay(k,e) (dashed lines) calculated with CPT for different 
values of the el-ph coupling A, using r\ = 0.05t. All other parameters as indicated in the figures. The truncation errors are 
A < (a) 9.1 x 10 -5 , (b) 9.0 x 10~ 5 , (c) 1.2 x 1CT 7 , (d) 2.4 x 10" 6 . The vertical lines correspond to variational diagonalization 
results for the polaron and bipolaron band dispersions (0) + -E T (fe) (dashed) and E^(k) (solid), respectively^ 



near the adiabatic regime due to the larger mass of the 
oscillators. 

In principle, the spectrum also contains coherent ex- 
cited states which are separated from the lowest-energy 
band by less than the phonon energy to. However, owing 
to the rather complex structure of the spectrum in the 
two-electron case, they are difficult to distinguish from 
the other contributions. A direct calculation of excited 
states in the Holstein model with one electron has re- 
cently been presented by Barisic^li Finally, the relation 
between A± and is very similar to Zu = 4. 

As we increase the el-ph coupling, the bipolaron dis- 
persion collapses to an extremely narrow band [Fig.^b)]. 
This cross over is again associated with a significant loss 
of spectral weight. At k = 0, for example, we find a 
reduction from about 0.50 at A = 0.5 to about 0.14 at 
A = 0.75. Increasing A further to 1.25, we finally arrive 
at a bipolaron band with W'/W ~ 10 -4 and a spectral 
weight of less than 0.03 at k = 0. Similar to Fig. ^d), 
the spectrum displays several bands equally spaced by us, 
which belong to states with one or more phonons excited. 
Moreover, the polaron and bipolaron bands are well sepa- 
rated, and the incoherent contributions dominate at large 



k. 



The agreement between the bipolaron band dispersion 
and E^(k) in Fig. is again very good. Similar to 
ZJ = 4, the condition for a small bipolaron is given by 
A > 0.5, so that CPT yields very accurate results. In 
contrast, the two-polaron energy E*(0) + E*(k) fits less 
well to the corresponding bands in the spectral function. 
We attribute this difference to the antiadiabatic regime 
(Fig. to the stronger retardation effects for ZU = 1. 
As a consequence, the polaron state is more extended 
below the small-polaron cross over occurring at A = 1 
(see, e.g., I), leading to a stronger residual interaction on 
a finite cluster, which also manifests itself in the CPT 
results. In contrast, for ZJ — 4, the lattice distortions 
around the electrons are very localized, and the two po- 
larons are almost independent. Above A = 1, i.e., in 
the small-polaron regime, the two-polaron dispersion for 
ZU = 1 again follows closely a two-polaron-like feature in 
the spectrum [Figs.EJc) and (d)]. 
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FIG. 4: (color online) Spectral functions Af(k,e) (solid lines) and Ay(k,e) (dashed lines) calculated with CPT for different 
values of the el-ph coupling A, using r\ — 0.05t. All other parameters as indicated in the figures. The truncation errors are 
A < (a) 3.3 x 10~ 5 , (b) 2.0 x 1CT 5 , (c) 9.9 x 10~ 6 , (d) 6.7 x f(T 7 . The vertical lines correspond to variational diagonalization 
results for the polaron and bipolaron band dispersions (0) + -E T (fe) (dashed) and E^(k) (solid), respectively^ 



2. U = 4 

So far, we have only presented results for U = 0, for 
which a bipolaron state is always favored. However, in 
materials such as the cuprates or the manganites, strong 
local correlations hinder the carriers from forming onsite 
bipolarons even for strong el-ph coupling. To model such 
effects, we therefore consider here a finite value of the 
el-el repulsion U = 4. 

In the case of two electrons with opposite spin, the 
Lanczos results for the cluster Green function converge 
faster as a function of iV p h for U > compared to U = 
as a result of the reduced effective el-el interaction. This 
is fortunate, since it allows us to use slightly larger clus- 
ters, thereby partly compensating for the increased finite- 
size effects due to the spatially more extended ground 
state in the weak-coupling regime. 

From the general discussion in Sec. |H] we expect the 
ground state to consist of two weakly bound polarons for 
2Ep < U, and a cross over to a bipolaron state at a crit- 
ical value of the el-ph interaction A. In the antiadiabatic 
limit, the latter is determined by 2Ep — U (i.e., A = 1 
for the case considered here) for weak coupling, and by 



4Ep = U for strong coupling. 

In Fig. we present the results for the spectral func- 
tion, again for A = 0.5-1.25. For weak coupling A = 0.5 
[Fig. Ufa)], the most striking difference to the U = 
case discussed above is the fact that there appears only 
one band at low energies. Together with the incoher- 
ent contributions, and taking into account the doubling 
of the number of carriers leading to a shift of energies, 
the spectrum bears a close resemblance to that of a sin- 
gle polaron with the same parameters (Fig. 3 of I). This 
is also underlined by the polaron and bipolaron band 
dispersions shown in Fig. Ufa), which are almost identi- 
cal throughout the Brillouin zone, and lie just below the 
corresponding band in the spectral function. In partic- 
ular, the band displays the typical flattening at large k, 
where the low-energy excitations have mostly phononic 
character. Furthermore, owing to the finite onsite repul- 
sion, the low-energy band in A± is very similar to that in 
Aij since, for finite U and weak el-ph coupling, the sin- 
glet ground state consists of two weakly bound polarons. 
Consequently, the singlet and triplet state have compara- 
ble energies, although the spectral weight in Ay is again 
very small near k = 0. 
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FIG. 5: Comparison of the spectral function A^(Q, e) calculated with ED and CPT, respectively, for different clusters sizes N, 
using n = 0.05t. The crosses correspond to the VDM result for the bipolaron energy £^(0)4- 



For A = 0.75 [Fig.EJb)], the ground state of the system 
is still given by two polarons, and the spectrum is almost 
indistinguishable from A = 0.5. In the present case, the 
condition for the existence of an intersite bipolaron is 
expected to lie between the weak- and strong coupling 
results U < 2Ep and U < 4EpH However, owing to 
its small binding energy, the intersite state is difficult 
to distinguish from the two-polaron state in the spectral 
function. 

At A = 1.0 [Fig. Efc)], the band in A T begins to split. 
Although the energy difference between the polaron and 
bipolaron band dispersions is still relatively small near 
k = 0, an excitation gap clearly emerges at larger k. 
Finally, at A = 1.25, two distinct bands with similar 
spectral weight have formed, which agree very well with 
E^(k) and (0) + E^(k), respectively. Interestingly, 
the band in the triplet spectral function A± lies notice- 
ably higher than the polaron band in Aj. Thus, for the 
parameters considered here, two polarons of opposite spin 
can lower their total energy by occupying the same lat- 
tice site, which is just the mechanism behind bipolaron 
formation. 

The abovementioned discrepancies between the bipo- 
laron band dispersion E^(k) obtained by Shawish and 
the band in Af are a result of finite-size effects in the CPT 
calculations. The latter become smaller with increasing 
coupling A together with the size of the bipolaron, and for 
A = 1.25 we find a very good agreement [Fig. EJd)]. To 
illustrate this point, we compare in Fig. [3] the spectral 
function Af at k = 0, A = 0.5, and for different clus- 
ter sizes N, calculated using ED with periodic bound- 
ary conditions (left column) and CPT (right column), 
respectively. The results reveal that for weak coupling 
and intermediate U, ED is superior to CPT concerning 
the convergence of the peak positions with respect to 



system size. This is not surprising as CPT is based on 
a strong-coupling expansion in the hopping termi Here, 
the el-el and el-ph interaction are both of about the same 
magnitude as the hopping, so that the method does not 
work as well as for U = 0. 

For U > 0, finite-size effects in both, CPT and ED, are 
larger due to the extended bipolaron state which exists 
for weak coupling. Similar to the one-electron case dis- 
cussed in I, deviations from the exact results due to the 
finite cluster size are usually smallest for k = 0, while 
they become larger with increasing k. Although in Fig.0 
the positions of the peaks in the CPT spectral function 
are slightly less accurate than in the case of ED, the 
weights of the excited states resemble more closely the 
results in the thermodynamic limit. 

Finally, for U > 4, the cross over to a small bipolaron 
occurs at even larger values of A. Apart from the change 
of the critical coupling, the physics is not altered signif- 
icantly. Therefore, we have restricted our discussion of 
the spectral function to U < 4, but some results for the 
bipolaron band dispersion at U = 8 will be presented 
below. 



3. Bipolaron band dispersion 

The bipolaron band dispersion E'^(k) has been cal- 
culated before by Wellein et alM and WeiBe et aZi 32 i 41 
for small phonon frequencies uJ = 0.4 and ZJ = 0.5, 
respectively. Remarkably, for parameters U > and 
A > such that the effective interaction U e g = 
[Eq. (J3J], they found a renormalized, free-particle dis- 
persion relation^Siii In this section, we wish to extend 
these considerations to the case 57 = 1, and to infinite 
systems. While the narrowing effect due to el-ph interac- 
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FIG. 6: Bipolaron dispersion E'^(k) as a function of the 
wavevector kw— Also shown is the bare tight-binding disper- 
sion for nearest-neighbor hopping, and a fit to the results for 
U — 4, A = 0.5 using a dispersion for nearest- and next- 
nearest neighbor hopping (see text). All curves have been 
scaled to the interval [—1,0], with the actual bandwidths given 
in the legend. 



tion has been discussed above, here we focus on the form 
of the band. 

Owing to the limited energy resolution and finite-size 
effects in the CPT results shown above, we use the 
more accurate data from the VDM. In Fig. we show 
Shawish'si^ results for the bipolaron energy as a func- 
tion of k, for different values of U and A. To permit a 
direct comparison, we have scaled all curves to the in- 
terval [—1,0], with the actual bandwidths given in the 
legend. 

We begin with the regime of a strongly bound small 
bipolaron. To this end we consider the case U = 
and A = 1.25. The corresponding band resembles quite 
closely a cosine dispersion, with some deviations being 
visible around k = tt/2. A different behavior is found 
for finite U = 4, as well as weak coupling A = 0.5. For 
these parameters, which favor a ground state with two 
polarons [see Fig.0{c)], the form of the band is remark- 
ably different from a simple tight-binding dispersion for 
nearest- neighbor hopping. This is still true for A = 1, 
although a trend towards a cosine dispersion is visible. 
For even larger U — 8, the noncosine-like form persists 
even for A = 1. 

It is worth mentioning the great similarity of the re- 
sults for U = 4 and U = 8 in the weak-coupling regime, 
which follows from the fact that once the small (onsite) 
bipolaron state is energetically unfavorable for the two 
electrons due to the Coulomb repulsion, a further in- 
crease of the latter has very little effect. On top of that, 
the intersite bipolaron state which exists for U > has 
a very small binding energy, so that the band dispersion 
is almost the same as that of two polarons. 

To identify the origin of the deviations from a free- 
electron band, we also included in Fig. El a fit of a free- 
electron model with nearest and next-nearest neighbor 



hopping to the band for U = 4 and A = 0.5, which 
yields an amplitude t' w 0.6t for two-site hopping pro- 
cesses. As proposed by Wellein et al.fi^ the importance 
of long-range hopping for the band dispersion of a single 
polaron may be due to a residual polaron-phonon inter- 
action, with the phonons and the polaron residing on 
different sites. Since we find substantial deviations of 
the bipolaron band from a cosine dispersion only in the 
regime of two weakly bound polarons, it stands to reason 
to assume the same underlying mechanism. 

Finally, we would like to comment on the fact that de- 
spite U c g — for [7 = 4 and A = 1 in Fig. we do not 
have a simple cosine band, in contrast to the findings of 
Wellein et al£& and Weifie et aL^L which have have been 
attributed to the formation of an intersite bipolaron^i In 
contrast, here we observe noncosine-like behavior even in 
the regime where an intersite state exists. These differ- 
ences are expected to be a result of the larger value of 
the phonon fre quen cy (here uJ = 1, while aJ = 0.4 and 0.5 
in Refs.l3fl andl4ll respectively), leading to a noticeable 
reduction of retardation effects. Moreover, since the crit- 
ical coupling A c decreases as 57 — > 0, t he bip olaron is more 
strongly bound in the work of Refs. 30 41J, thereby sup- 
pressing the abovementioned nonlocal phonon-polaron 
interaction. Further work along these lines is highly de- 
sirable to understand the dependence of the bipolaron 
band dispersion on the phonon frequency in the regime 
uJ < 1. 



V. CONCLUSIONS 

We have presented a detailed study of the one-electron 
spectral function of the Holstein-Hubbard model with 
two electrons of either the same or opposite spin. The 
method employed here is cluster perturbation theory to- 
gether with the Lanczos method, which represents a ver- 
satile and fast approach. 

As a function of the electron-phonon and electron- 
electron interaction strength, polaron and bipolaron 
states manifest themselves as quasiparticle bands, and 
results have been compared to accurate data for the bipo- 
laron energy dispersion. For weak coupling and/or in- 
termediate to strong Hubbard repulsion, finite-size ef- 
fects are visible, but are much smaller than in previous 
work restricted to small clusters. The major advantage 
of the present method is that the spectrum can be ob- 
tained at any point in fc-space, even when using clusters 
with only a few lattice sites for which enough phonon 
states can be kept in the calculation. This has allowed 
us to investigate, for the first time, the dispersion and the 
spectral weight of the quasiparticle features throughout 
the Brillouin zone. The results and their dependence on 
the model parameters have been discussed, and a perfect 
agreement has been found with the physical picture of 
the Holstein-Hubbard bipolaron emerging from previous 
work. A comparison of the bipolaron dispersion with a 
simple tight-binding band has revealed an important con- 
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tribution from next-nearest-neighbor hopping processes 
in the regime of a weakly bound state. 

Finally, the adiabatic regime of small phonon frequen- 
cies, characteristic of many real materials, remains an 
interesting and demanding open issue for future work. 
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